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TECHNICAL PAPER 


A PRELIMINARY STUDY OF NUMERICAL SIMULATION OF THERMOSOLUTAL 
CONVECTION OF INTEREST TO CRYSTAL GROWTH 

I. INTRODUCTION 


Thermosolutal (or “double-diffusive”) convection shall refer, in this work, to buoyancy-induced 
convection where there are two or more density components (e.g., heat and solute). The most interesting 
cases are when the diffusivities of the separate components are very different from each other. In the most 
common example (heat and salt), the diffusivities differ by a factor of about 100, heat having the higher 
diffusivity. This results in some perhaps surprising phenomena — in particular, the instability of a quiescent 
fluid which is bottom-heavy and which has no horizontal density gradients. 

The two simplest situations (continuing with the heat/solute example) are those in which there are 
no horizontal gradients of either component, and where one component contributes to a stable stratification 
while the other component contributes to an unstable stratification. For the purposes of discussion, it is 
assumed that the overall density increases downward (i.e., the fluid is bottom-heavy). If no double-diffusive 
effects are considered, a fluid element displaced vertically would have a force exerted upon it which would 
restore it to its initial position. Due to dissipation, it would arrive there with less momentum (and 
buoyancy) than it started with. Thus, vertical motions are damped and no instability would arise. However, 
if the more slowly-diffusing component contributes to a top-heavy state, a phenomenon, called “salt 
fingers,” can occur. Since the diffusivity of heat is very large compared with that of the solute, the vertically 
displaced element tends to “feel” the destabilizing buoyancy effects of the solute rather than the effects of 
temperature. That is, the fluid element finds that it is quickly heated (or cooled) to the approximate tem- 
perature of its new environment, but that it has a different density than its environment due to its solutal 
concentration, which remains approximately the same as ihe original value. Thus, its vertical displacement 
is accelerated rather than damped. This instability consists of vertically long convection cells with small 
horizontal extent (hence the name “salt fingers”). The second situation is when the fluid is stable with 
respect to the more slowly diffusing component (solute), and unstable with respect to temperature (i.e., 
hot at the bottom). A vertically displaced fluid element will be buoyantly restored toward its initial position 
(due to the stable stratification of the solute) and will arrive there with essentially the same solutal con- 
centration with which it left; but because of the relatively high diffusivity of heat, after it arrives back at 
its original position the fluid element finds itself either warmer (if the initial displacement was downward) 
or cooler (if the initial displacement was upward) than it was originally. Hence, an amplifying oscillation 
develops, most often called the “oscillatory” or “diffusive” instability. The appearance of this type of 
convection is usually that of horizontal layers of convection cells; but this form is much more complex 
than the finger convection, and the number of layers within a given vertical depth depends upon initial and 
boundary conditions and may be time dependent [ 1 ] . 

It is reasonable to expect that thermosolutal effects may be important in crystal growing situations 
where there exist density gradients in the fluid feed material due to both temperature and solutal concen- 
tration gradients, and that in some cases any theoretical model which neglects these effects is omitting some 
fundamental physics. The crystal growing situation that shall be used as our “typical” system here is the 
Bridgman-Stockbarger method, whereby a substance in a long cylindrical ampoule is held vertical and moved 
downward through a furnace with highly controlled temperature conditions. The substance is thus solid 
(crystalline) below a hot melt. A basic axial (vertical) temperature that is buoyantly stable is thus imposed 
upon the melt, the lower surface being approximately isothermal (at the melting point temperature). The 



substance under process is often made of two components. One of the components (the “solute”) tends 
to be rejected at the interface and incorporated into the solid less readily than the other substance (the 
“solvent”). Thus, an axial gradient in solute may exist in the system. If the solute is lighter than the solvent, 
it is speculated that the “finger” process may give rise to that simplest form of thermosolutal convection. 
A linear analysis of this situation is performed in Reference 2. Additionally, radial (horizontal) gradients in 
both temperature and concentration may induce convection, which may be expected to be influenced by 
the double-diffusive nature of the system. For example, even when the solute is heavier than the solvent, 
horizontal temperature gradients combined with the stable concentration profile may give rise to horizontal 
layers such as those observed by Thorpe et al. [3] and Chen et al. [41 in a heat/salt system. The presence 
of these layers, with very narrow regions of large gradients, would make the numerical modeling of the 
fluid motions quite challenging [5,6,7] . The modeling of the fluid motions is important in order to predict 
the characteristics of the crystal itself - for example, the concentration of solute (“dopant”) in the crystal 
may have important radial variations if cellular convection is present, and axial variations if the flow is time- 
dependent. 

Rather than give an extensive review of research in thermosolutal convection in general, this report 
refers the reader to the references herein, especially Chen and Johnson [8] , which is a review article of a 
conference on the subject in the spring of 1983. Most of the work has been done under the auspices of 
oceanographic research beginning in the 1960’s (actually, the first work on it was Stommelet al. , [ 9 ] ), and a 
rather sparse scattering of experimental, observational, analytical, and numerical work has been done. 
Some of the work that has been done on numerical modeling of thermosolutal convection, with particular 
emphasis placed upon applicability to crystal growth, is reviewed here. 

The linear stability of a quiescent state has usually been handled by purely analytical methods, the 
computer being used only for the purpose of displaying the results. Although not discussed here, some of 
this work is included in the bibliography [10,11,12]. Veronis [13,14] and Sani [15] were the first to 
compute fully nonlinear thermosolutal convection. They did so for the oscillatory situation between hori- 
zontal flat plates (the infinite channel) which has limited applicability to crystal growth. However, it is 
interesting to note that in some instances finite amplitude disturbances were found to be unstable even 
when the thermal driving was subcritical — i.e., when the system was stable to the infinitesimal disturbances 
of linear theory. Huppert and Moore [16] extended the work of Veronis to give more details about the 
type of flow expected as a function of the parameter combinations. Their work and the work of Moore 
et al. [17] indicate the complexity of the system, even when values of the diffusivities involved are not as 
extremely different as would be expected in physical situations. The finger system is less complex, but still 
difficult to model numerically for all but the highly idealized systems. Straus [18] used a Galerkin tech- 
nique to study fully nonlinear finger convection between horizontal plates. He made progress toward 
reconciling the fact that observations of the cells indicated much shorter wavelength than that predicted 
by linear theory. He found that only small-scale motions are themselves stable, and that the most stable 
convection cell has a wavelength which is approximately the same as that which maximizes the flux of 
solute. The present author has been unsuccessful at finding previous calculations of either type of convec- 
tion in an enclosed container. 

As indicated earlier, there have been previous calculations performed of a system with horizontal 
temperature gradients imposed upon a fluid in a container which is initially stably stratified with a salt 
gradient [5,6,7]. These calculations met with limited success. From their figures it can be seen that the 
spatial resolution was not adequate, since a spatial computational mode (“wiggles”) is apparent. The calcu- 
lations were successful, however, at predicting the horizontal layering that occurred in the laboratory experi- 
ments, including the correct wavelength. The difficulty lies in the fact that the horizontal convection cells 
have very strong shear between them, with accompanying strong gradients over narrow regions in the tem- 
perature and concentration fields. Thus, high spatial resolution is required to resolve all the fields (especially 
concentration) — resolution which is not practical to use, except perhaps on today’s super-computers (not 
available to the previous workers). 
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Calculations of crystal growth fluid dynamics have not yet been made with full allowances for 
thermosolutal effects, although there is work presently being done in this area by groups at MIT (under 
the direction of R. A. Brown) and at the National Bureau of Standards. (F. M. Carlson [19] at Clarkson 
College has considered purely thermal convection.) At NBS, G. McFadden and S. Coriell (Schaefer et al., 
[20] ) have performed calculations of melt flow including thermosolutal effects where there is a crystal 
“growing,” but the presence of sidewalls and interface curvature are neglected (i.e., a horizontally infinite 
channel is considered). A finite difference method is used, and a computational domain is chosen such that 
one complete wavelength in the horizontal is resolved, and such that the upper horizontal boundary is far 
j enough from the flow (which is present due to the unstable boundary layer at the interface - the lower 

| boundary) that there is approximately no influence from it. A growth rate is imposed, and due to a nonzero 

; segregation coefficient, (light) solute is rejected at the interface. The heat flow due to conduction is calcu- 
lated in the solid crystal below. This work is an extension of the infinite channel work noted above, the 
primary difference being that the basic temperature and solute fields are exponential, rather than linear, due 
to the imposed growth rate and boundary conditions. “Salt finger” convection is predicted for certain 

i parameter combinations, and the NBS group has investigated such questions as the vertical extent of the 

convection, the solute and heat fluxes, and the horizontal wavelength at which the convection cells split 
into two. No attempt has been made to determine the wavelength at which the convection is most vigorous 
in some sense (such as maximum solute transport). The recent published work by the MIT group [21,22] 
has made the transition from calculations of crystal growth including convection due to temperature 
j gradients only to that with solutal buoyancy effects as well. Calculations are performed in a cylindrical 
domain, and include the heat flow in the solid crystal as well as the calculation of the interface shape. The 
sidewall boundary conditions on temperature in the fluid domain are those of no heat flux for a given 
: distance near the interface, and isothermal above that. The convection studied is due to the radial 

(horizontal) temperature gradients that are present because of this discontinuous boundary condition; the 
convection is either enhanced or suppressed by the solute, depending on whether it is light or heavy. This 
\ is not thermosolutal convection in the sense described above, since there is no qualitative change in the flow 
due to the presence of the solute. The fact that the flow (except for the parallel flow due to the crystal 
growing) is due entirely to the discontinuous boundary condition, makes detailed study of it of questionable 
relevance to actual crystal growth, where heat flow in the ampoule may cause the boundary conditions in 
! the fluid region to induce an entirely different (surely quantitatively if not qualitatively) flow than that 
studied by the MIT group. It should also be noted that the diffusion parameters and density variations used 
by the MIT group are not nearly as extreme as those of the actual physical situation, and thus thermosolutal 
effects such as fingers or horizontal layering were effectively excluded. 

The present work is considered to be preliminary results of an effort to study thermosolutal convec- 
tion of particular interest to crystal growth by means of numerical models. Finite difference methods are 
used, the models being adapted from those used to study laboratory experiments and the infinite channel in 
geophysical fluid dynamics [23,24,25] . Since explicit time differencing is used mainly, computer resources 
required are rather large — the primary reason these results are “preliminary.” Section 4 (Finger Convection 
in the Infinite Channel) describes some progress toward making the codes more efficient. After the improve- 
ments are completed, research will be oriented toward areas that results have shown to be the most 
promising for further study. There are three areas studied: (1) horizontal heating of a stably stratified 

solution; (2) finger convection in a box, in some cases with the inclusion of a horizontal temperature 
gradient on the bottom surface; and (3) finger convection in infinite channel geometry. 


II. SIDE-HEATED CASES 


The calculations reported in this section relate to the laboratory results of Thorpe et al. [3] , Chen 
et al. [4] , Wirtz et al. [5] , and Hart [26] — and, to a lesser extent, to Ruddick and Turner [27] . In the first 
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four works, a box is filled with a salt solution, the concentration of which decreases linearly upward (i.e., 
the fluid is made bottom-heavy). One of the side walls is then heated, and the other held at a fixed (lower) 
temperature. In References 3 and 26, the heating is done very slowly; one purpose being to allow a basic 
flow to become established with a horizontal gradient in solute concentration, which approximately com- 
pensates for the horizontal gradient in temperature. In the vertically infinite situation, there exists a known 
steady basic state — an exact solution to the full equations where there are no vertical variations in the 
velocities, temperature, or the solute gradient. This situation can be approximated by a very tall, narrow slot 
in the laboratory. The experiments by Chen et al. and Wirtz et al. (these two works are by the same group) 
had the temperature of the heated side wall rise very fast (with an e-folding time of about 3 min). The 
problem in this case is that of the stability of a time-dependent flow near a single vertical wall. The work by 
Ruddick and Turner [27] was a set of laboratory experiments where a sugar solution was beside a salt 
solution (both stably stratified), and the partition between them was carefully removed. In all of the side- 
heated cases, appropriate parameter combinations resulted in a flow with horizontal layers, where each 
layer consisted of a convection cell of thermally direct rotation, with very strong shear between them. The 
development of cells began near the heated side wall, and the cells would eventually extend across the width 
of the volume. The theory of Hart [26,28] explained very well the slowly-heated tall slot, and a simple 
theory explains the vertical depth of the layers in the rapidly side-heated case [4] . Chen [29] numerically 
studied the linear stability of the time-dependent parallel flow in the latter case, with good success in pre- 
dicting the critical wavelength and Rayleigh number. Even though the flow is essentially two-dimensional, 
numerical modeling of the fully nonlinear flow is difficult (even when cases with only a few layers are 
considered), because the strong shear between layers and the low diffusivity of the solute results in the need 
for very high resolution. Another difficulty that can arise is due to the observations (and predictions by the 
theory of Holyer [30]) that the layers themselves may become unstable to finger and diffusive convection. 
The finger convection, especially, places an even more severe constraint on the resolution required to 
describe the flow. An attempt to avoid cases where this might occur was made in the present work, although 
it may be seen that the resolution, while somewhat greater than that used in previous works (except, in a 
sense, for Wertz et al., [5] ), was evidently inadequate in some cases, and perhaps with finer resolution these 
detailed features would have been predicted. 

This report shall consider two-dimensional flow of a fluid in which the Boussinesq approximation is 
valid, z is taken to be the vertical coordinate and y to be the horizontal coordinate. If the two components 
of the momentum equation are cross-differentiated, a vorticity equation is derived: 

-3iMf d*ar 1 r2 , m 

3t 3z dy dy dz p Q dy 

where v is the kinematic viscosity and where the stream function is related to the vorticity by: 


V 2 i> = f 


( 2 ) 


The velocity components are related to the stream function by: 
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3z 
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The equations describing the evolution of the temperature and solute concentration fields are, respectively: 
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(5) 


where k is the thermal diffusivity and D is the solutal diffusivity. Density shall be taken to be a linear 
function of temperature and solute concentration: 


p = p o [l-a(T-T o ) + 0(C-C o )] 


( 6 ) 


where a and 0 are the thermal and solutal coefficients of expansion, respectively. In this section, boundary 
conditions on temperature shall be taken to be of the Dirichlet type on the side walls (i.e., they are assumed 
to be perfect conductors), and the horizontal end walls will be perfectly insulating. In some cases the tem- 
perature of the hot wall will be brought up exponentially (with time), and in others it will be fixed through- 
out the integration. The side walls, and in some cases the end walls as well, shall be taken to be zero-flux 
walls with respect to the solute. There are some cases where the end walls have fixed solutal concentration 
imposed upon them. No-slip boundary conditions are used in all cases. 

If the equations are non-dimensionalized, the dimensionless parameters that appear are the thermal 
and solutal Rayleigh numbers, the Prandtl number, and the Schmidt number. Aspect ratio also appears 
after the spatial domain is defined. These are defined: 


g « ATL 

Rcl'-r 

1 VK 
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(9) 

Sc = v/D 

(10) 

height 
y~ ... 

01) 


weight 


5 



A finite difference technique is used to solve the equations. The domain (a vertical plane) is 
sectioned into a grid - not necessarily of uniform grid element size, but one where the resolution may 
gradually increase near the boundaries. The equations are finite differenced in explicit centered form, 
although because of the grid stretching the spatial differencing is not truly centered (i.e., second order). 
(A variable transformation scheme is not used.) Details of the technique are in Miller and Gall [23,24] . The 
technique was shown to be valid by reproducing in detail a laboratory transition curve between axisymmetric 
and non-axisymmetric flow for the baroclinic annulus. Furthermore, axisymmetric calculations of 
thermally-driven flow (rotating and non-rotating) in an annulus agree quite well with those of previous 
workers. The necessary addition to the Miller and Gall code was simply to add the solute concentration field 
to the problem. After calculating the temperature, concentration, and (interior) vorticity fields at a given 
time step, the finite differenced equation (2) is solved either by direct Gaussian elimination (for the 
abruptly-heated cases) or by successive line over-relaxation (SLOR). (The latter technique proved to be too 
inefficient when the stream function had rapid change over a small region.) The no-slip boundary condition 
on the stream function (normal derivative equals zero) is imposed when the boundary value for vorticity is 
then calculated. The initial conditions are that there is no motion in the fluid volume. Computer time 
required on the VAX 780 for the resolutions here is about 30 min per 1000 time steps. 

The rapidly-heated side wall case shall be discussed first. This case uses the no-flux boundary condi- 
tion on solutal concentration on the entire boundary. (Note that when this is the case, the only true steady 
state is one where there are no solutal gradients in the fluid volume — eventually, the flow will mix the 
solute throughout the fluid.) Chen et al. [4] reports that the initial rise of temperature of the side wall 
could be modeled by an exponential function AT(1 - e -t / T ), where r is about 2 min. Thus, the code was 
modified to allow for this; i.e., the temperature that was imposed upon the hot wall (the other one being 
fixed for all time) rose exponentially, with an exponential rise time of either 2 min (cases 1 and 2) or 3 min 
(case 3). The initial condition for the interior was of constant temperature — that of the “cold” wall. 

Chen et al. points out that the distance L that a fluid particle, within a basic density gradient 
0 O = |3 3C Q /3z, receiving a temperature increase of AT, will rise before it is in a neutral buoyancy environ- 
ment, is: 


a AT 

L = . 

$0 

It is expected that the vertical spacing between layers will be of this scale, and this is verified experimentally. 
Chen et al. uses this L as the length scale in defining the dimensionless parameters, and so will we. The 
dimensions of the computational domain in this report are referred to in terms of this L. 

It is noted that the experimental results were for many layers (on the order of 100). It is not realistic 
to numerically model this situation under the present circumstances (explicit methods, small computer). 
Thus, the attempt was to model cases where one would expect five or fewer layers. Results are given here 
for cases with one to four layers. Note that Wirtz et al. [5] performed limited calculations (the integrations 
were performed only until the layers were just beginning to develop) for a case with about eight layers, and 
Wirtz [7] performed calculations for cases with only one and two layers. The relevant numerical and 
physical parameters used here are given in Table 1. 

The vertical extent of case 1 (Fig. 1) was 5L, and the horizontal extent was 3L. (The dimensions 
were 3 cm by 5 cm for this case, and 2 cm by 3 cm for all other cases in this section.) The diffusivities were 
the actual ones of heat and salt, and the viscosity was that of water. The resultant Prandtl number was about 

7, Sc = 714, Ra T = 5 x 10 4 , and Ra s = 1.5 x 10 5 . Initially, as the temperature difference is brought up 
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TABLE 1 . RELEVANT NUMERICAL AND PHYSICAL PARAMETERS 


Case 


Ra s 

Height 

Pr 

Sc 

Time step 

Grid 

Total Time 

1 

50,000 

150,000 

5L* 

7.14 

714 

0.03 sec 

32x31 

285 sec 

2 

10,000 

30,000 

3L 

7.14 

714 

0.03 sec 

32x31 

450 sec 

3 

40,000 

40,000 

1L 

7.14 

714 

0.03 sec 

32x31 

600 sec 

4 

21,000 

63,000 

3L 

7 

10 

0.05 sec 

30x28 

300 sec 

5 

4,200 

12,600 

3L 

7 

625 

0.05 sec 

30x28 

2700 sec 

6 

2,940 

8,820 

3L 

1 

10 

0.05 sec 

30x28 

215 sec 


*Height is given in terms of the buoyancy length L = a AT /0 Q . Dimensional size of L was 1 cm for all 
cases, and viscosity was 1 centistoke. 


from zero there is a single cell throughout the volume (Fig. 1 a), although it may be seen in the figure that 
cells are beginning to break off in the lower and upper comers near the heated wall. Note that the time has 

been so short (35 sec) and the velocities so small (about 10'^ cm-sec'*) that a fluid particle has moved only 
a very small fraction of the total height. At a later time (142 sec) an intense cell has developed in the lower, 
hot corner. This is as observed in the experiments. A weaker cell is evident near the upper boundary. 
After 285 sec, two layers between the corner circulations are beginning to form, resulting in four layers. 
Although five layers were expected, it is noted that the experimental result was that the thickness of the 
layers was only approximately L, and also that the comer circulations change the basic temperature and 
solutal fields from the original one. At this point in the integrations, however, it is evident from the solutal 
concentration field (and even more so for the field of vertical density derivative, not shown) that resolution 
is not adequate, because there are wiggles which correspond to every other grid point. An upstream differ- 
encing scheme was tried for the solutal equation, and although this resulted in smoother fields, it did not 
predict the layering and thus was even more undesirable. 

Case 2 (Fig. 2) was for thermal and solutal Rayleigh numbers one-fifth those in case 1 . (This is a 
sub-critical case, according to the laboratory results.) There is a strong circulation in the lower warm corner, 
but layers that would be identifiable in the lab are not apparent. Yet still there are indications in the solutal 
concentration field that resolution is inadequate. Case 3 (Fig. 3) is for a case with height equal to 1L, and 

the thermal Rayleigh number is equal to 4 x 10^. The circulation is qualitatively like one computed by 
Wirtz [7] , although it is integrated here for a longer time period. Obviously, the resolution is inadequate, 
although the temperature field and the general shape of the flow are quite reasonable. 

Calculations were also performed with more idealized parameters, looking for interesting flows that 
could be well resolved using the present techniques and computer. Case 4 (Fig. 4) has Ra-p = 21,000 and 

height 3L. Prandtl and Schmidt numbers were 7 and 10, respectively. Besides the much higher solute 
diffusivity, this case also differed from the previous ones in that the temperature of the hot and cold walls 
were fixed for all time, and that the initial condition for temperature was that of a linear gradient between 
the two walls. Thus, the problem was symmetric with respect to the hot and cold wall. The boundary 
condition on solute concentration was again no-flux for all surfaces, and the flow shown is quasi-steady. 
There are only two layers present, one circulation centered near the lower hot wall and another near the 
upper cold wall. The circulation is such that the concentration and temperature fields are adjusted so that 
the density contours in the interior are nearly horizontal. Note that qualitatively this flow is very similar 
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STREAMFUNCTION 

STREAMFUNCTION 

STREAMFUNCTION 

MAXIMUM - 

1.36555E— 03 

MAXIMUM * 

7.62230E— 03 

MAXIMUM - 

9.38488E— 03 

MINIMUM 

- 1.46924E-08 

MINIMUM 

— 6.18865E-- 04 

MINIMUM 

—1 .1 1702E—03 

INCREMENT * 

1. 50000 E— 04 

INCREMENT - 

1.00000E— 03 

INCREMENT = 

1 . 50000 E— 03 


TEMPERATURE 
MAXIMUM - 2.2158 

MINIMUM = — 1.90191E-04 

INCREMENT = 0.30000 


CONCENTRATION 
MAXIMUM * 9.84237E— 04 

MINIMUM = — 1.03723E— 03 
INCREMENT * 3.00000E-04 


Figure 1 . Results of calculations for case 1 of Section II (see Table 1 ). (a) Stream function at 35 sec 
(b) Stream function at 142 sec. (c), (d), and (e): Stream function, temperature, and concentration 
at 285 sec. Stream function value should be divided by 500 in Figures 1 through 12. 






STREAM? UNCTION CONCENTRATION TEMPERATURE 

MAXIMUM « 5.76839E— 03 MAXIMUM = 1 96 726 E -04 MAXIMUM = 0 46502 

MINIMUM * — 2 .47488 E— 04 MINIMUM = -2.05857E-04 MINIMUM = 7 70147E-05 

INCREMENT = 8.00000E-04 INCREMENT = 5.00000E-05 INCREMENT = 6.00000E-02 



Figure 2. Results of calculations for case 2 of Section II (Table 1) at 450 sec. 


STREAM? UNCTION 
MAXIMUM - 7.7337 

MINIMUM = -0.16325 
INCREMENT = 1.0000 


CONCENTRATION 
MAXIMUM * 6.321 26E— 05 

MINIMUM = -9.69310E— 05 
INCREMENT = 2.00000E-05 


TEMPERATURE 
MAXIMUM = 0 68866 

MINIMUM = - 1. 44892 E -03 
INCREMENT = 8.00000E-02 



Figure 3. Results of calculations for case 3 of Section II (Table 1) at 600 sec. 
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to one of the cases computed by Wirtz [7] , the main difference being that Wirtz used a Schmidt number of 
1000 (that for salt), and he had a resolution problem because of it. Some of his comments regarding the 
physical significance of the fine structure of the solutal field are not well substantiated, since the fine 
structure is not resolved by the grid and may be numerical rather than physical. 

Case 5 (Figs. 5 and 6) is for a Schmidt number of 625 and Raj = 4200, other parameters being the 

same as for case 4. The boundary conditions on solute at the top and bottom are of fixed concentration. 
(This means that a steady state, which retains the solutal stratification, is possible.) At 300 sec there are 
again the two circulations at the bottom near the warm wall and at the top near the cool wall. Here, these 
circulations are much more separated (i.e., of smaller vertical extent), a consequence of the lower solutal 
diffusivity. Since solute is more nearly conserved, the buoyant restoring force is effectively stronger, and a 
fluid particle rising (sinking) because it is heated (cooled) does not go as far before turning inward. Note 
that again the interior isopycnals (Fig. 5d) are approximately horizontal. 

After a very long time of integration (2700 sec), the solution has not completely converged, but 
changes occurring are quantitative only and are occurring very slowly. There are three convective layers 
(Fig. 6), the middle one becoming substantial in the contour plot at about 1200 sec. It is hypothesized 
here, although not argued on precise theoretical grounds (this would be future work), that the 3-layer flow 
resulted out of the instability of the flow pictured in Figure 5. The solutal and thermal Rayleigh numbers, 
as defined by Hart [28] for this case, are very near the critical curve (on the unstable side) for the vertically 
infinite slot. It is reasonable to hypothesize the possibility that a basic state would be able to develop 
when one is so near the critical curve; and that although the present system is very far from being an infinite 
slot, the critical curve for this system may not be too far from the idealized one. Note that the fields are 
fairly well resolved, although the solute concentration field does exhibit a weak computational mode, 
indicating that resolution should be better. 

Perhaps the most interesting case is depicted in Figures 7 and 8. This is case 6, and is for Raj = 

2940, Pr = 1, 8c = 10, and height equal to 3L. The concentration is fixed at the top and bottom, as in 
case 5. The flow is vacillating between a single-cell thermally direct (counter-clockwise) circulation and a 
thermally indirect circulation. These circulations and the transitions between them are shown in the 
sequence of stream function plots in Figure 7. Note that a positive stream function value implies a counte- 
clockwise flow and that a negative stream function value indicates a clockwise flow. The transition between 
single cell flow is characterized by a break-down of the single cell into multiple cells and an expansion and 
eventual merging of cells of the opposite rotation (opposite from the previous single cell flow). The 
thermally direct cells form and expand from the lower hot and upper cold corners. The thermally indirect 
cells do likewise from the interior. The pattern shown here repeated itself several times, and the amplitude 
of the oscillation (in, say, kinetic energy with time) seemed to be rather constant. Further work is needed 
for a definitive demonstration that this is indeed vacillation (i.e., that it will never become a steady flow). 

Figure 8 shows the temperature and concentration fields at time 203 sec and density field at time 
203 and 200 (similar to that at 212) sec. The temperature contours remain rather vertical; the concentra- 
tion field approximately compensates for the horizontal temperature gradient so that the net horizontal 
density gradient is small. It is evidently the small deviation from no horizontal gradient, however, that 
drives the vacillation. Note the opposite tilt of the isopycnals in Figure 8c and 8d. In Figure 8d, the cir- 
culation is near its counter-clockwise peak, and the isopycnals are nearly horizontal. The circulation advec- 
tion quickly causes a tilt of the interior isopycnals so that there is a reversing torque (Fig. 8c), and an 
opposite circulation develops which eventually tilts the isopycnals in the opposite direction (not shown), 
again reversing the torque. This process is perpetuated by the boundary conditions, which cause the per- 
sistance of a horizontal density gradient in the lower hot and upper cold comers. It may be hypothesized 
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Figure 6. Results of calculations for case 5 of Section II (Table 1) at 2700 sec. 

that there exists, as a solution to the equations, a weak therjnally direct circulation, but that it is not seen 
because it is unstable. It is also possible, as previously suggested, that if this integration were carried out 
much further the vacillation would go away. This is a specific area where further work would be interesting. 

In conclusion, although some very interesting and informative work can be done with the present 
spatial resolution for idealized parameters, to accurately model the laboratory experiments will require 
much higher resolution. This work indicates that if similar types of flow exists in crystal growing systems 
(and it is reasonable to expect so), then the resolution required to model the flow in the melt and to 
therefore accurately calculate such quantities as dopant segregation will be quite large; it is likely that 
super-computers such as a Cray or the Cyber 205 will prove to be useful (in fact, necessary) for this task. 


111. FINGER CONVECTION IN A BOX 


The same model, described in the previous section, was used to calculate some flows where the 
vertical gradient in the temperature field was in the stable direction and the solute gradient was destabilizing. 
The boundary conditions on temperature and concentration were of fixed value on the upper and lower 
boundaries and of no flux on the side boundaries. No-slip boundaries are assumed, and the same dimensions 
(2 cm by 3 cm) and viscosity as in the previous section were used. When there were no horizontal gradients 
imposed, a very small random perturbation in the concentration and temperature fields was added to the 
initial condition, which was that of diffusive equilibrium and no flow. Some of the cases have horizontal 
temperature gradients on the bottom, which is analogous in a certain sense to a curved interface in crystal 
growth. An interface may be curved, for example, when the thermal conductivity of the solid below the 
melt is significantly less than that of the melt. The heat must flow downward, and in order to do so some 
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Figure 7. Calculated stream function for case 6 of Section II (Table 1) at various times, (a) 206 sec, (b) 209 sec 
(c) 2 1 0 sec, (d) 2 1 1 sec, (e) 2 1 3 .5 sec, (0 2 1 6 sec, (g) 2 1 7 sec, and (h) 2 1 8 sec. 
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Figure 8. (a), (b), and (c): Calculated temperature, concentration, and density fields for case 
of Section II at 203 sec. (d) Calculated density field at 200 sec. 






of it must be conducted through the ampoule (side wall). This requires the presence of horizontal tem- 
perature gradients and hence a curved interface, since it is approximately isothermal. 

Although no theory has been developed for the stability of the finger situation when side walls are 
present, the critical solutal Rayleigh number can be approximated by consulting the infinite channel theory. 

Case 1 (Fig. 9a, b, c) is for Rag = -1.2 x 10^, Raj = -8.3 x 10^, Pr = 1, and Sc = 10, where the total height 
(3 cm) was used for the length scale in the Rayleigh numbers. Using the free-slip infinite channel theory 

(Baines and Gill, [31]), the critical solutal Rayleigh number, given the other parameters, is -8.3 x 10^. 
This case is expected to exhibit an instability, but it would not be very strong. The flow depicted was 
integrated for 1050 sec, and is still evolving (growing in amplitude). Being relatively near the critical curve 
resulted in the growth rate of the instability being small and hence taking a long time for the convection to 
become the amplitude seen in the figure. There are one and one-half complete horizontal wavelengths 
contained within the volume. The advection effects upon the solute are strong, but the effects upon the 
temperature field are rather small. 

Case 2 (Fig. 9e, f, g) was for the same parameters as the previous case, except that the Schmidt 
number was 25. Hence, the critical Rag was -3.3 x 10^, and the Rag was about four times critical. The 

flow was much stronger, the wavelength smaller, and the instability arrived faster. The fields shown in the 
figure for this case are also still evolving. 

If there is a horizontal gradient imposed upon the lower surface and no solutal effects are considered, 
a single thermally direct cell is induced (case 3, Fig. 10). Because the horizontal driving is only on the lower 
surface and because of the imposed vertical gradient, the circulation is confined to the lower region of the 
volume. Some calculations were performed to test the effect that the addition of this circulation would 
have upon the pure salt finger cases described above. Case 4 (Fig. 1 1) is with the horizontal gradient added 
to case 1, and case 5 (Fig. 12) is with the same horizontal gradient added to case 2. The addition of this 
lower circulation seemed to “kick off’ the finger convection faster than with the quiescent basic state. 
The overturning, which tends to stabilize the basic state in the pure thermal convection case, did not seem 
to perform an analogous function here. None of the cases here (except the pure thermal case) were inte- 
grated to steady state, and it is possible that cases 1 and 2 would eventually have become more vigorous 
than cases 4 and 5. It is important to point out that the horizontal temperature gradient was imposed by 
lowering the temperature of one of the lower comers, rather than by increasing the temperature of the 
opposite corner. The amount was by about 10 percent of the total vertical temperature difference, and 
hence if the opposite had been done, the magnitude of the thermal Rayleigh number would have been 
effectively decreased by about 1 0 percent on the warmed side — which would have increased the instability 
by itself, even without the presence of the induced circulation. We emphasize that such was not the case 
here. 


The computer resources, required for the simulations here, were rather large (similar to those of the 
previous section) since they were performed before improvements in the computer code, such as implicit 
differencing of the diffusion terms, were made. Some of these improvements were made on the infinite 
channel code, and the improvement in efficiency (reviewed in the next section) was quite dramatic. 

It would be quite useful to have detailed experimental results of the configuration studied in this 
section, for the purpose of rigorous validation of computer codes. Previous laboratory studies of salt fingers 
have been made by producing the fingers at an interface between hot salty water above cold fresh water, 
rather than having linear gradients of heat and salt through the entire depth. To construct a basic state of 
the latter type would be much more difficult. 
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Figure 9. (a), (b), and (c): Results of calculations for case 1 of Section III at 1050 sec. (See text for details.) 
(d), (e), and (f): Results of calculations for case 1 of Section III at 700 sec. 
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Figure 12. Results of calculations for case 5 of Section III at 175 sec 





In conclusion, there is evidence here that the always-present horizontal temperature gradients in the 
Bridgman-Stockbarger system may actually contribute to, rather than suppress, finger-type convection. 
Thus, any full simulation of the crystal growth situation where a light solute becomes concentrated at the 
interface must allow for the possibility of this relatively fine-scale convection. (That is, the resolution must 
be rather fine.) It is noted also that, due to the observation of the experimental works relating to the 
previous section, fingers may appear between horizontal convective layers when the overall solutal gradient 
is stable, and that this may occur in a crystal growing situation as well. 


IV. FINGER CONVECTION IN THE INFINITE CHANNEL 


The calculations presented here were performed as part of an effort to monitor a contract to demon- 
strate the performance of a spectral model on areas of relevance to crystal growth. The performance of the 
finite difference model here will eventually be compared with the performance of the spectral model. Free- 
slip conditions were taken on the horizontal boundaries, because the spectral model was based upon this 
condition and converting the spectral model to no-slip would not have been a trivial effort. The problem is 
essentially the same as that considered by Straus [18] , although Straus assumed very small solutal diffusion 
compared with thermal and viscous diffusion, and thus simplified the equations. We have made no such 
simplification. The model here was based upon that of Miller [25] which is similar to the one used in the 
previous section with the following exceptions: boundary conditions on the vertical end walls are taken to 
be periodic. Solution of the Poisson equation [equation (2)] was by Fourier analysis in the horizontal and 
direct Gaussian elimination in the vertical. Improvements upon the model described in Miller [25] are as 
follows. The finite differencing of the vertical part of the temperature field was done implicitly (Crank- 
Nicholson), greatly improving the efficiency of the model for the PR < 1 cases reported here. (For the 
Pr = 0.1 cases, the allowable time step was increased by a factor of three.) Further improvement would 
have been achieved if the viscous diffusion term had also been differenced implicitly, and if ADI methods 
had been implemented to allow implicit differencing in the horizontal. Allowance was made for the time 
steps to be different for the vorticity, temperature, and concentration prediction equations. In the steady 
state situations studied here, this resulted in considerably fewer time steps required for convergence. This 
is because the small diffusivity of solute would otherwise require integration for a very long time in order 
to achieve equilibration. The improvement factor in computing the solutal Nusselt number to the degree 
of accuracy in Table 2 is about equal to the Schmidt number (10). These improvements were suggested by 
previous work of Dr. Glyn Roberts of Roberts Associates, Inc. (Roberts et al., [32] ). The number of time 
steps for the coarser resolution reported here was about 2000 for “complete” convergence (constant solutal 
Nusselt number to 5 digits), which took about 4 min of CPU time on the UNIVAC 1180. Before the 
improvements were made to the code, this calculation would have taken about 2 hr. 

Initial conditions were of diffusive equilibrium, but with a perturbation added to the solutal con- 
centration field. The perturbation was in the form of a complete sine wave in the horizontal and a half sine 
wave in the vertical, of size 0.01 . 

Results are shown for Pr = 0.1 and Sc = 10. The length scale used in defining the Rayleigh numbers 
in these cases is the vertical depth of the channel. When the solutal Rayleigh number is super-critical enough 
that the flow is substantially nonlinear, this results in a moderately difficult problem for the numerical 
models, with a factor of 100 separating the two extreme diffusivities - although the diffusivities are not as 
extreme as in a true crystal growing situation. In that case, the Prandtl number would be about 0.01 (for 
liquid metals), and the Schmidt number would be between 100 and 1000. The assumption of periodic end 
wall conditions is also a great simplification over a more realistic domain (such as that in the previous 
section); only one horizontal wavelength is modeled here and hence fewer grid points are required. The 
results are reported by plotting the fields, plotting energy versus time, and by giving the solutal Nusselt 
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TABLE 2. NUSSELT NUMBER 


Case 

Raj 

Ra s 

P 

Sc 

L 1 

Grid 

Nu s 2 

la 

-1000 

- 90.9 

0.1 

10. 

2.828 

12 x 10 

3.548 

lb 

-1000 

- 90.9 

0.1 

10. 

2.828 

16x24 

3.544 

2 

-1000 

-100 

0.1 

10. 

2.828 

12 x 10 

3.68 

3 

-1000 

-100 

0.1 

10. 

1.95 

12 x 10 

3.74 (3.74) 

4 

-1000 

-200 

0.1 

10. 

1.58 

12 x 10 

4.83 (4.85) 

5 

-1000 

-300 

0.1 

10. 

1.47 

12 x 10 

5.58 (5.62) 

6 

-1000 

-400 

0.1 

10. 

1.40 

12 x 10 

6.17 


1 . L is horizontal length of the computational domain, in terms of the vertical depth. 

2. Numbers in parentheses are results from Straus [ 18] . 


number, which is the ratio of the resultant solutal flux through the boundaries to that of pure diffusion. 
Note that calculations have been performed only for cases where the resultant flow is steady. If the solutal 
Rayleigh number is increased enough, periodicity and chaos can ensue (S. A. Orszag, private communication, 
1983). 


The first case shown here is one that is only moderately super-critical. The thermal Rayleigh number 
was -1000 and the solutal Rayleigh number was -90.9. The critical Rag for this situation was about -76. 

The horizontal wavelength was 8*'^, which is the critical wavelength at the critical Rag. Figure 13 shows 

the vertical resolution of the two cases, along with the resultant average concentration profiles. The hori- 
zontal resolution was 12 intervals for the first case (Fig. 13a) and 16 intervals for the other. The number of 
vertical intervals was 10 and 24, respectively (the size of which was smoothly varied within the vertical 
extent of the domain). No qualitative differences are apparent between the two resolutions, and quantita- 
tive differences are small. In particular, the solutal Nusselt number (which would be expected to be rather 
sensitive to the resolution differences) is 3.548 and 3.544, respectively. Thus, the seemingly coarse resolu- 
tion of 1 2 x 10 performs rather well, although the reader is reminded that the grid is stretched in the vertical, 
so that the resolution near the boundary is increased from that of a constant grid. In Figure 14 are plots of 
the fields at steady state for the coarser resolution case. The flow consists of convection cell extending the 
entire height of the domain. The temperature field deviates only slightly from the conduction state, and 
the solutal field is highly deformed. Note from Figure 13 that the horizontally averaged concentration field 
actually has a reversed vertically gradient in the middle of the domain. 

Figure 15 shows the contour plots (at the coarse resolution) for the case with Rag = 300, which 

results in a more highly nonlinear flow. Note that the solutal perturbation is much larger, and that the 
convection cells are less symmetric with respect to the vertical. From the contour plot of concentration, it 
is evident that the resolution is beginning to be suspect, but the Nuc comparison with Straus [18] is still 
quite good. ° 

In summary, great improvements in computer usage efficiency can be made with the use of implicit 
methods in the diftusion terms and, for steady state solutions, with the use of different time steps for the 
dilferent prediction equations. Furthermore, calculations of this configuration demonstrate the accuracy 





CONCENTRATION CONCENTRATION 

Figure 13. Calculated horizontally-averaged solutal concentration field as a function of height of two 
resolutions of case 1 of Section IV (Table 2). The discrete points indicate the vertical 

grid used in the calculations. 

of the grid stretching method and difference scheme, since the results compare very well with the spectral 
calculations of Straus [18], even for the low resolution grid. As previously noted, however, these calcula- 
tions are for a single convection cell, and for multiple cells the resolution would have to be increased 
appropriately. 


V. CONCLUSIONS AND DISCUSSION 


The calculations here indicate both the potential and shortcomings of finite difference methods for 
modeling thermosolutal convection. The shortcomings are not unique to finite difference methods, but 
would exist for finite element and spectral techniques as well — although the degree of shortcomings of each 
method is a subject for debate. The difficult aspect of modeling the full nonlinear flow is that thermo- 
solutal convection tends to be of fine scale, due to a low-diffusivity component. This results in the need for 
many degrees of freedom in order to accurately describe the flow. Furthermore, the existence of higher- 
diffusivity components results in time step restrictions (which may be quite severe, due to the fine resolu- 
tion), unless implicit methods are used. (Note that all three methods almost invariably use finite differencing 
in time.) Implicit methods in the finite difference context have the advantage of being easier to implement 
than for the other two methods. 

The flows calculated here certainly reinforce the expectation that double-diffusive convection is 
important in crystal growing systems, although much more work is needed before this is shown to be a fact. 
Certainly it is reasonable to expect that thermosolutal convection is important when there is a tendency for 
a light component to be present near the bottom of the liquid, even when the liquid is cooled from below. 
This is the “salt finger” situation. Furthermore, even when the solutal contribution is to make the system 
bottom-heavy, horizontal temperature gradients, if they are strong enough, can conceivably drive a flow 
which causes horizontal layering and subsequent overturning of the solute (c.f., Section II). Experimental 
results have indicated that finger convection can then occur within the individual convection cells. Compu- 
tations of this phenomenon were not performed, due to the very high resolution that would be required. 


23 






O CM 

z 2 

3 g 

s°' 

<< 

LLt 

cc 


GO 


ii 


cc < 


LU 

cc 

D 

k r 
< o 
cc II 

8 ! < 

2 



2 < 









Figure 15. Contour plots of the calculated fields for case 5 of Section IV (Table 2). 









The main aspects of the crystal growing situations that were neglected here are: (1) possibly curved 
lower boundary; (2) a basic growth flow through the system, which results in exponential rather than linear 
background states; (3) cylindrical geometric effects; and (4) three-dimensional effects. To fully include all 
of these effects would indeed be a giant step forward in the modelling of crystal growth. The first three 
effects are included in the crystal growth modelling efforts of the groups mentioned in the introduction. 
(Including fully three-dimensional effects may not be a realistic goal, even with the use of the latest genera- 
tion of super-computers — at least when thermosolutal convection is present.) The primary common 
weakness of those efforts mentioned is that they did not fully include the thermosolutal effects studied 
here. It is hoped that the present work will assist crystal growth modellers by pointing out the extent of 
some of the resolution and other numerical problems facing them in order to have a complete description 
of the fluid dynamics of the melt, and hence of the crystal itself. 
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